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ABSTRACT 

i 

For Courant numbers larger than one and cell Reynolds num- 
bers larger than two, oscillations and in some cases instabili- 
ties are typically found with implicit numerical solutions of 
the fluid dynamics equations. This behavior has sometimes been 
associated with the loss of diagonal dominance of the coeffi- 
cient matrix. It is shown here that these problems can in fact 
be related to the choice of the spatial differences, with the 
resulting instability related to aliasing or non-linear interac- 
tion. Appropriate "filtering" can reduce the intensity of these 
oscillations and in some cases possibly eliminate the ’ instability . 
These filtering procedures are equivalent to a weighted average 
of conservation and non-conservation differencing. The entire 
spectrum of filtered equations retains a three-point character 
as well as second-order spatial accuracy. Burgers equation has 
been considered as a model. Several filters are examined in de- 
tail, and smooth solutions have been obtained for extremely large 
cell Reynolds numbers . , 
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1. Introduction 


Three-point finite-difference discretization has generally 
been used to approximate the partial differential equations of 
fluid mechanics. Explicit formulations are typically restricted 
by the "linear" stability conditions 

c = uAt/Ax^l and R = uAx/v^2 , (1) 

where c is the courant number; R c is' the cell Reynolds number; 
u is a reference velocity; At, Ax are the temporal and spatial 
increments, respectively. 

Implicit central difference formulations designated by 

{ ICD ) are linearly unconditionally stable; however, the diagonal 

dominance of the tridiagonal inversion matrix is assured only 

if conditions (1) are satisfied, if conditions (1) are violated, 

i.e., there is a loss of diagonal dominance, error growth is 

possible and spurious oscillations are observed. Recent- 

ly a difference approximation, denoted by (KR) , that insures 

linear stability and still maintains diagonal dominance, even 

( 2 ) 

when conditions (1) have been violated has been proposed. ' For 

♦ 

linear systems, solutions have been obtained for c»l and R c »2. 
Unfortunately, it has been shown that the KR diagonally dominant 
formulation also exhibits instabilities for the non-linear Bur- 
gers equation, when conditions (1) are violated; moreover, for 
the same conditions, stable solutions are sometimes obtained with 
the non-diagonally dominant central differences.^^ Therefore, 
it may be concluded that non-linear instability (aliasing) and not 
diagonal dominance is the primary reason for error growth when 
conditions (1) are violated . This may explain the difficulties 
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t 


for c»l that have been encountered with implicit Navier-Stokes 
finite-difference solvers. since, lack of diagonal dominance 
leads to spurious oscillations^, for non-linear systems this 
enhances the aliasing effect. 

, , (4) 

Previous studies on non-linear instability by Phillips, 
Arakawa,^ Piacsek and Williams ^ and others have been pre- 
sented for explicit schemes; these, satisfy the CPL condition C,<1. 
However, the lack of conservation of certain quadratic quantities 
over the whole domain still leads to an aliasing error growth. 

A smoothing procedure, or alternatively, proper spatial differ- 
encing (filtering) can help to eliminate this instability, in 
the past, first-order accurate upwind differencing has been used 
extensively, as it has sufficient numerical viscosity to suppress 
these instabilities. This corresponds to a severe filter. Other 
schemes, most often used for hyperbolic systems, e.g., the MacCor- 
mack scheme and variants of the Lax-Wendroff method, do not ap- 
pear to encounter this problem. This is, perhaps, due to the 
presence of 'increased amounts of numerical viscosity, either in- 
herent or added artificially to the finite .difference system. 

In this paper, the Burgers equation is re-evaluated for both 
the I CD and KR techniques. It is shown that (a) for the linear’ 
Burgers , equation, with the procedure of reference (2), there is no 
instability even if conditions (1) are grossly violated. This 
confirms the linear analysis of reference (2) and the earlier re- 
sults of reference (3) . Error growth, with the ICD technique, 
due to the loss of diagonal dominance, is possible, but it is shown 
that this can be avoided in most cases for the linear Burgers 
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equation? (b) for the non-linear Burgers equation, it is shown 
that instabilities can appear with the ICD and KR formulations, 
even if the conditions (1) are not simultaneously violated. 

This would again indicate that diagonal dominance is not respon- 
sible for the resulting error growth? (c) for the non-linear 
Burgers equation, if smoothing or filtering of the non-linear 
term is applied with the method of reference (2) , there is ap- 
parently no instability even if conditions' (1) are violated? in 
addition, spurious oscillations are considerably diminished. 

The magnitude of the oscillations associated with the lack of 
diagonal dominance is closely related to the resulting non-linear 
instability. When these high frequency oscillations are filtered, 
even maintaining second-order accuracy , the aliasing effect is re- 
duced. With a filtered ICD formulation, aliasing can also be con- 
trolled? however, the shock may convect to the boundary of the 
domain. This is true regardless of the boundary location. It 
will be shown that this solution with the shock located at the 


boundary is in fact an alternate solution of the ICD difference 


equations. This non-uniqueness is due to the lack of .conservation 
N 

of I. xx. over the domain, and occurs only when the cell Reynolds 

j=l 3 A 

number R = . 

c v 

Finally, the linear as well as non-linear stability of the 


ICD difference equations is evaluated around a given initial state* 
This analysis confirms that filtering procedures reduce the magni- 


tude of oscillations and have a positive influence on the non-linear 
instability. Smoothing procedures of the type utilized by Shuman ^ 
and investigated in more detail by Shapiro ^ are not considered in 


the present paper. ORIGINAL PAGE IS 

’ ‘ OF POOR QUALITY 
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2. Analysis 

Burgers equation describing a fixed shock wave is given as: 
u t + (u-*s)u x = vu xx . ‘ (2) 

The initial and boundary conditions are taken to be: 
u(x,0) = 1 for -5sx<0 

u (0, 0) « .5 . (3a) 

u(x,0) = 0 for 0<x£5 


and 

u(x,t) = 1 for x ~-5 

(3b) 

u(x,t) = 0 for x=5 

The exact solution is given to a good approximation in the steady state 
hy * 

u » | [l-tanh 

All the solutions discussed in the present paper are obtained with 
51 equally spaced grid points? Ax = 0.2. Convergence was assumed 
when differences in the values of u between 100 iterations were 
such that |u k+100 -u*!<10~ 6 ; k denotes the iteration number. Many 
calculations were also run with more severe convergence conditions. 

The results were unchanged. 

Numerous calculations of the linear version of equation (2), 
where the coefficient of the convective term is treated as a con- 

\ 

stant, were obtained with both the implicit central differencing 
(IQD) and the Khosla-Rubin (KR) scheme. The restrictions on cell 

i.i 

Reynolds number and Courant number, necessary for diagonal domi- 
nance, were grossly violated. In every case a converged solution, 
for both schemes, was obtained. For the same conditions, with the 
non-linear Burgers equation, 10 iterations were performed at each 
time step. For- the KR scheme the calculations diverged. The solu- 
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tions approach a final converged solution to within 10~ and then 
exhibit a rapid divergence. For cell Reynolds numbers greater 
than two the solutions obtained with c<l or for central differ- 
ences in some cases with c>l, exhibit, the expected oscillatory 
behavior. As the steady state is approached, aliasing effects 
begin to accumulate, and in most cases .for c>l, R c >2, and even for 
some where c<l, R c >2, the solutions diverge for large times. In 
some of these cases, where R c is slightly larger than 2, a solu- 
tion with central differencing can sometimes be obtained, (e.g., 
v=l/24, or R c =2.4)? the ICD scheme has a somewhat larger transient 
numerical viscosity than does the KR formulation. The above des- 
cription of linear convergence and non-linear divergence, with 
equivalent values of c and R q , indicates that it is probably not 
the loss of diagonal dominance but the aliasing effect, which is 
enhanced by spurious oscillations, that is the reason for error 
growth . 

In the following sections, we will show that the source of the 
aliasing error lies in the form of spatial discretization. For 
conventional non-conservation form, non-linear instability arises 
with cell Reynolds numbers greater than two? this result is essen- 
tially independent of the value of the Courant number. For explicit 

schemes, this non-linear instability was first observed by Phil- 
(4) 

Ifps. One of the remedies he recommended to suppress such error 

growth was the use of smoothing of the short wave length components. 
Shuman devised a smoothing operator that is quite often used in 
■numerical weather prediction calculations. Alternatively, modified 
spatial differencing can suppress the high frequency modes. The fil- 
tering procedure is applied • here . Three different filters, including 
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one due to Shuman will be investigate'd in conjunction with the ICD 
and KR schemes. The aliasing error growth can then be eliminated 
so that a steady converged solution is obtained. First-order 
temporal and second-order spatial accuracy of both the implicit 
schemes under investigation are retained with all filters. 

2.1 Smoothing and Filtering 

In the present section, we will discuss a series of filters 
that can be used to suppress the aliasing error growth. Only 
three-point second-order accurate filters are considered. These 
are given by 



u . , , +u . . +ku . 
1+1 i-l i 

2+k 


(4) 


It can be shown that these filters (4) correspond to a weighted 
average of conservation and non-conservation difference equations. 
The filters (4) are second-order accurate, with a truncation error 
Ax u /(2+k) . Several specific cases of interest are: 

XX 

i) k = 0 is conservation form or a trapezoidal filter; 

ii) k = ® is non-conservation form (no filter) ; 

iii) k = 2 is the Shuman filter; 

iv) k - -1 is a special filter to be discussed later. 

In general, k can take on any value greater than negative two. How 
ever, the optimum choice is governed by the minimum amount of numer 
ical viscosity required to suppress the aliasing error growth. The 
V^lue k = -1 is significant, as will be shown in a later section; 
for large flow Reynolds number with k = -1 the effective cell Rey- 


nolds number is always less than or equal to two. 

It must be emphasized that there is no general way of arriving 
at an optimum value of k for any given problem. This can only be 
obtained by numerical experimentation and some physical insight in-. , . 
to the flow characteristics, in addition, the degree of filter- 
ing will depend upon the choice of the mesh. For example, for a 
flat plate boundary layer governed by the Blasius equation, the 
non-conservation form of the equations does not lead to stable 
solutions for grid spacings larger than 2 fsee Appendix]. Also, 
for non-uniform grids having large grid spacings in the outer por- 
tion of the boundary layer the solutions exhibit oscillations and 
are, in general, rather poor. The use of filtering eliminates 
many of these problems. With a trapezoidal filter, converged 
solutions can be obtained with a minimum of grid points within the 
boundary layer. For non-uniform grids the oscillations in the 
regions of large mesh size are also reduced or eliminated? however, 
the second-order numerical viscosity introduced by the filter af- 
fects the accuracy of the surface shear stress. This can be re- 
duced by taking a smaller grid near the wall. For moderate uniform 
grids the artificial viscosity always leads to a less accurate shear 
stress as compared w-ith the non-conservation solution. Therefore 

for boundary layers, it would appear that different filters (k val- 
ues) should be used near and far from the surface. 

Once again, it should be stressed that these filters are not 
applied to smooth the solution after a given time step, but are used 
to recast the nonlinear coefficient of the u term in equation (2) . 

It is in this respect that the present investigation is different 
from the work of Shuman^, and similar to that of the Arakawa^ 
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and William ^ schemes, in what follows# we will examine certain 
of these filters for two finite-difference schemes, i , e, , ICD and 


KR. 


( 1 ) 


2.1.1 Central Differencing (ICD) 

Implicit central differencing as is well-known'"' is uncon- 
ditionally stable. The general form of the difference equations 
•with the filter (4) is as follows: 

At 


u n+1 *= u n - 
U j j 2 Ax 


r/ u n +i + 

P+l 3-1 

+ku? +1 > 

Li 

l 2-t-k 

y 

1 2 


,,,n+l n+1* 

(u j+i - U j-1 } 


+ ^ < u £i - ^j 1 + *T- 1 i ) 


(5) 


In earlier studies, the smoothing character of the conservation 
difference equations (k = 0) has been demonstrated by many inves- 
tigators (1.3) and will not be repeated here. The resulting oscil- 
lations or instability for R c >2 with non-conservation form 
(k - co or no filter) are also well-documented. 

Numerical solutions of (5) have been obtained for k = «, 0, 

2, see Table 1. A typical example is given by the conditions 
At = 6.0 and v = 1/96. This corresponds to a Courant number of 
15 and a cell Reynolds number of 9.6. The non-conservation equa- 
tions do not lead to any solution. Although a converged solution 
is obtained for k = 0,2, the shock wave is found to convect to the 
upstream boundary for these conditions. It will be shown that the 
finite difference equations ((5)with k « o) , in fact, has two solu- 
tions and that one is consistent with this convective behavior. 

Since the filter (4) is simply a linear combination of conservation 
and non-conservation equations, we will investigate the possible non 
uniqueness for these forms of the governing equation (5) . These are 
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and 


du. 

3t“ 


+ u 3+i~ u i_ J-. , . i (u , _ 2u+u , 

dt ^ 2 2 R ^ u j+l j j-l' 

(conservation) (6a) 

+ (U^Js) (Uj+i-Uj.!) = | ( U . +l -2u j +u.. 1 ) , 

(nonr-conservation) (6b) 

Where t is redefined as t/2Ax and R = The governing equation 

. V 

• if CO - 

(2) satisfies the conservation law j— ( Judx) = 0. The discrete 

~ CO 

forms of Burgers equation (6a) and (6b) should also satisfy this 
Conservation property. Summing over all grid points, we find 


d 

at 


N 

(E u.) 

j=l 3 


U 2 +U N-l r ‘ i 


^■ U 2~ U Nr>l r R^ 


and 


(conservation) 


d_ 

at 


N 

(S u.) 

j=l ? 


(I,. 

' 2 


r) ^ 1 ~ u 2“ U N-1^ 


(nonr-conservation) 

It can be seen that unless the symmetry condition is strictly en- 

j=l, the possibility of a "spurious 11 unsymmetric 
solution, for which .u 2 ^u N _ 1 r-2/R = 0, exists for the conservation 
form (6a). The solution for which the wave convects to the boundary, 
in fact, satisfies this latter ’condition. It has been found, on the 
other hand, that with a small At (Courant number <1) the symmetric 
solution is usually recovered. In the next section, it will be 
shown that for R£2 the symmetric solution is unique. For R>2, the 
non-symmetric solution also exists and is also stable. For non- 


forced, i.e, U 2 +U N 
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conservation form (6b) ; only the symmetric solution is possible. 


2.1.2 KR Scheme 


This scheme, as introduced in Reference 2, is diagonally 
dominant and unconditionally stable for all R and retains second- 
order spatial and temporal accuracy of the convective derivative. 
It is given by 


and 




+ ^ D ? +1 
Ax^ 3 


(7 a) 


for 


n+1 n At , n+1 . « .n+1 n+L At r f^i+l ^ n n 

_u - " (u ^ 2 ’ > + TKZ -H) Dn 


+ D^ 1 

Ax^ 


j+1 j 


2 Ax 


(7b) 


for u j <% 

where dJ = uS.,- 2u.+ u^ , . 

3 D+l 3 3-1 

For the KR formulation, the filters are only used to modify the 
, r non-linear coefficient of as shown in (7) ? the implicit portion 
of the convective term is always given by an upwind differencing 
form. The resulting equations are second-order accurate in A;t. 

Unlike the ICD results, for v = ^ and At = 6.0, both the 
trapezoidal and Shuman filters do lead to converged solutions des- 
cribing a thin symmetric shock. The results are shown in table 1. 
Although the convergence condition is satisfied, there is some indir- 
cation of a creeping motion of the shock wave. After 3000 time 
steps this movement was still less than the prescribed tolerance. 

* • I 

The infinite time behavior was not considered. For explicit schemes, 
this problem is virtually undetectable, as the Courant numbers are 
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limited by the CFL stability condition. The allowable At is much 
smaller than that considered here and therefore examination of 
the largo time history would require an extraordinary number of 
time steps . 

2 ,2 Non-linear and Linear Stability 

In the present section, we examine the stability of the 

I CD finite difference equivalent of equation (2) about a given 

initial state. Since the source of aliasing error lies in the 

*■ 

spatial discretization; the problem will be examined for a semi- 
discrete system. The underlying idea is that if the spatial dif- 
ferencing leads to a temporal amplification, then there should 
not be any time discretization that can eliminate this instabil- 
ity. For example, it will be shown, though not rigorously, that 
for non-conservative form and central differencing, the solution 
of equation (2) grows exponentially if the cell Reynolds number 
is greater than two. The stability will be examined for central- 
differencing for both non-conservative and conservative forms., 
The effects of other types of filters will also be indicated. 
2.2.1 Non-Conservation Form 

With central differencing, equation (2) becomes 


du . . 

ar + (n r h) - (u j+i" u j-i } = r <Vi“ 2u j‘ fu j-i ) 


( 8 ) 


r is a non-dimensional time. As it is difficult to find a closed 
form solution cf equation (8) for arbitrary j, only a few simple 
cases will be examined for small numbers of grid points N. 

(i) N ~ 3 : Since there are only three points, and from the boundary 

conditions we require u^ = 1 and u^ - 0, then the equation for u 2 


becomes 


du, 

A 

dT 


“ (u 9 - h) (l-2u 0 ) 


( 9 ) 
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The solution of equation (9) is 


u. 


u, (0) 


-(§ -l)t 
e + 


h (i 


~(f -Dt 

-e ) 


( 10 ) 


where U2(0) is the initial value of u 2 (t). Clearly the solution 
diverges for R>2. However, if ( 0 ) = %, as in. the present case. 
Equation (10) leads to u 2 -h for all t. However, in this case, 
there is a cancellation of two growing terms, and therefore the 
numerical solution can still be amplified by roundoff errors. 

This is seen to be the case for forward time marching in (9) . If 
we use an explicit method, we obtain 


u” +1 =u£(l+At ~ Tp + Jk (| -1) ) . (11) 

Clearly any small error will grow if the coefficient of u 2 is 
greater than unity. Thus, stability is possible only if R<2. 

For an implicit scheme with the convective coefficient (Uj~%) 
treated explicitly, then (9) becomes. 


u 


n+1 


<1 -» 

2 (1+|^) 


1+ At 


n 


2 At u 2 
R 


1+! 


( 12 ) 


Once again, the Neuman stability condition requires R<2 for a 
stable solution. If, we treat ( u 2 ~^) implicitly, the solution 
converges for almost all At and R. Hpwever, such a scheme is 
inconsistent with the’ differential equation (9) . 

(ii) N = 5: In this case, u, = 1 and u,. = 0. 

X O 


9 1 

"dt + = R ( u 3“ 2u 2 + D ' 

- i 

du ■ 

dt" + (u 3 ~%) (u 4 -u 2 ).- - (u 4 -2u 3 +u 2 ) , 

“ - u 3 , (u 4“^ ) - = R (u 3~ 2u 4 ) 


( 13 ) 


dt 


A closed form solution of equation (13) is possible, if we assume 
the shock to be symmetric about u 3 =%. Thus, we get 


u 2 +u 4 = 1 and 


du„ , 3 

dtT + ^R U 2 = 2R “ ^ 


(14a) 


du 4 , ,2 , . 1 , 

dt + R U 4 ~ 2R 4 

integration of Equation (14a) gives 


(14b) 


u 2 = u 2 (0) e R ■ * R - 


-(# -h) t 


(15) 


Once again the solution grows with time, unless R<4. From the 
two cases considered, it appears that if the mid-point symmetry is 
not fixed the solution will diverge . for R>2 while a converged 
symmetric solution may be possible for 2<R<4 if the mid-point 
value u 3 =% is fixed. The non-conservative calculations of Refer- 
ence (3) confirm the validity of this inequality for the cases, 
considered. 

2.2.2 Conservation Form 

A similar procedure can be carried out for the conservative 
difference form of Burgers equation. 

(i) N = 3 ; In this case, the governing system reduces to 

dU 2 2 _ 1 
dr“ + R u 2 ~ R 

Integration gives 

u 2 (t) = h + (u 2 (0)-%) exp (~2 T /R) . (16a) 

S.ignif icantly , this system leads to a steady converged state with 
u 2 ( T ) = h for all R. 
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N = 5: 


The governing equations are 
du Uo(Uo-l) i 

ar + -4“ - I < V 2 u 2 +1 » 


du 

dT 


3 . , , W 1 

+ (u 4 -u 2 ) —2 


= R (» 4 -2“ 3 +u 2> ' 


du 


u 3 (u 3- X) 1 


R ( U 3 -2U 4^ 


dT 2 

Two steady-state solutions are possible, 
u 2 +u 4 =l and u 4 -u 2 =’ 2/R 

The symmetric solution characterized by u 2 +u 4 =l is stable for all R 
The non-symmetric solution is given by: • • 


u 3 = l+yi^-4(4/R 2 -|) , 

2 

l+2u +4/R 

u 2 = , (16b) 

l+2u^-4/R 
U 4 = 4 

0 

This solution exists only for R>2.5. In order to investigate the 
stability properties, we perturb about (16b) and look for the solu- 
tions of the type, exp (It). This leads to the following dispersion 
relation 

* = r ± /-s - (2 U °-D 2 / 4 

e R e - 

Since \ always has a negative real part, the solution (16b) , in 

the steady state, is stable for R&2.5. It should be pointed out 

that the solution with enforced symmetry (i.e., u^Jz) is stable for 

all cell Reynolds numbers. This has been numerically tested for 
-5 5 

v = 10 or R=10 . The resulting solution has oscillations but is 
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stable. These oscillations can be eliminated or reduced in 
amplitude by applying the filters of the type discussed pre- 
viously * 

2.2,3 Optimal Filtering 

We recall that the different filters are defined by 
„ u . , ,+u . ,+ku. 

U .= . JH, .J-ri. 1 . (17) 

u j 2+k ' '■ n 

for k = 0, we recover the trapezoidal filter or conservation 
form; for k = °=, non- conservation form is recovered; k = 2 cor- 
responds to the Shuman filter and k - -1 is a noteworthy case. 

The finite-difference form of Burgers equation (2) with the 
filter (17) is a weighted average of non-conservation and con- 
servation equations and is given as 

du. u., .+u. i+ku. , 

_1 + J . . Jj] (u j+1 -u j _ 1 )= | 

(18) 

We shall examine the case where N = 3 in order to obtain an opti- 

* 

mum k value for (18) . For N = 3, (18) becomes 


du, 

dT 


l+ku 0 , 

^ 2+k ^ ~ R ( 1-2 u 2^ 


so that 


U 9 (t) = h + [u o {0) -%] exp ( 


k 2, 
2+k " r' T 


2 ' ' 2 L 2 

For large R, k>0 has a destabilizing influence and k^O has a 

stabilizing effect . We recall that the filtering introduces 

to u' 


(19) 


an artificial viscosity 


, when -2<k<0 . For k<-2, the 


2 (2+k) 

sign of this viscosity changes and consequently the filter is no 
longer useful. From the previous analysis of non-conservation or 
conservation solutions, it is known that oscillatory behavior oc- 
curs when R>2. These oscillations can be eliminated by reducing 
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the grid size and therefore the local value of R. The filter (17) 
can also accomplish this without grid reduction. Por large 
R{»2) a minimum amount of filtering is required in order to ob- 
tain smooth solutions. The degree of filtering, as characterized 
by- the value of k, should be such that the changes in u^ are con- 
fined to only one grid point. This amounts to incorporating arti- 
ficial viscosity such that the cell Reynolds number based on the 
"effective viscosity" does not exceed 2; the optimum choice for k 
can be seen by comparing (19) with (16a) such that R eff = 2 or 

1 14-k 

R 2+k 

• . (9) 

This relation has also been obtained by Cheng and Shubin' ' from 
different considerations. Results for several filters are pre- 
sented in the following section. 

3 ■ Results 

Numerical solutions using various filters were obtained for a 

variety of Courant and cell Reynolds numbers. These do confirm 

* 

the stability analysis of the previous section. Por example, for 
large cell Reynolds number and Courant number, a converged solution 
for conservation form is obtained by enforcing the symmetry condi- 
tion. This solution exhibits oscillations, characteristic of large 
cell Reynolds number, flows. Without enforced symmetry, the shock 
may convect to the boundary. This corresponds to the second solu- 
tion discussed previously, in 'the cases considered here, symmetric 
solutions are usually found for Courant numbers less than one. in 
these cases a completely symmetric solution is not achieved. The 
shock continues to move with an extremely small velocity. However, 
within a prescribed tolerance, the solution can be considered to be 


16 


ORIGINAL PAGE IS 
OP POOR QUALTTXi 


converged. Calculations for a variety of cell Reynolds numbers 
ranging from 2.4 to 50,000 and various filters characterized 
by -l^ksl, were carried out by Taverna and Busch, ^ For a 
given cell Reynolds number, an optimum filter was defined by 
a hiimimum error condition. Velocity (u) profiles with and with- 
out filters are shown on figure (1) . The optimum values as a 
function of R (figure (2)) are also reproduced from reference 
(10) . It should be noted that with an appropriate filter, oscil- 
lations characteristic' of large cell Reynolds numbers can be 
eliminated; formally, second order accuracy is still maintained. 

For R- ro , we note that k optimum - -1. 

4 , Summary 

1. Diagonal dominance problems associated with TCD methods 
can be eliminated by the KR scheme. Calculations with a linear 
Burgers equation confirm the analysis of reference 2. 

2. Stability problems arising in calculations with locally 
large cell Reynolds numbers are found only for non-linear equations 
and are due to the form of spatial differencing of the convective 
terms. This instability and associated oscillations can be elimin- 
ated by appropriate filtering, 

3. A non-uniqueness of the conservation form of the differ- 
ence equations is described. The second (non-physical) solution is 
encountered numerically only for large Courant numbers and R>2. 

4. Finally the results of 1-3 have been confirmed by numerical 
experimentation as well as some approximate stability analysis. 

5. The present analysis has been confined to the stationary 
solution of Burgers equation, although some results for boundary 
layers are given in the Appendix. For more general equations, the 






nature of the optimum filter may vary from that obtained here. 
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Blasius Equation 

The flat plate boundary layer in .similarity variables is 
governed by 


u p + fu T) = 0 

f l) = u 


(non-conservation form) 


+ (fu) ^ -u = 0 


(conservation form) 


The finite-difference form using weighted averaging with h - AT] 
is given by 

u i-g- 2 V u i-i f jfi +kf i u k£ j +£ l-i u "”1— i— . es o 

h 2 + 2(l+k)h j+1 2 (l+k)h u 3-l 1+k 

(A.l) 

clearly k = 0 and k = » lead to the conservation and non-conserva- 
tion forms. In order to examine the effect of filtering, the 
truncation error of equation (ft.l) is investigated. Taylor series 
expansion about the jth grid point lead to 

2 2 2 iv 

< 1+ 3 u+k) >U T)11 * (f+ = ~ ^ fU + 2 ^ (A,2) 

The second-order accuracy of the numerical scheme is retained. The 
additional truncation error arising out of weighted averaging is 
shown as coefficients of the convective as well as diffusive terms. 
It may be seen that the convective and diffusive modification can 
be made small near the surface by taking a fine grid. However, near 
the edge of the boundary layer, where large h values can lead to a 
deterioration of the solution resulting in oscillations, the numeri- 
cal viscosity of the filter can be quite large and thus the oscilla- 
tions are suppressed. The effect of filtering is therefore to in- 
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corporate damping where it is needed the most. The filtering ef- 
fect is largest for -l<k<0 for these equations. In the present 
context <?f filtering, conservation form provides a significant 
amount of artificial viscosity, so that solutions with large mesh 
sizes are possible. For example, a converged solution with h = 6 
(the boundary layer thickness is about 3.5) can be obtained; non- 
conservation solutions are no longer possible when h>2. This 
artificial viscosity provides a thickening of the boundary-layer 
and consequently a reduction of the wall shear. Smaller grids 
near the wall surface are required to eliminate this accuracy 
problem. Optimally, it would appear that non-conservation form 
should be used near the surface and conservation form in the outer 
portion of the boundary layer. This corresponds to a variable 
filter. 
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TABLE X I SOLUTION OF BURGERS EQUATION 
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FI6. la FILTERED BURGERS EQUATION: R *2.4 





FIG. lb FILTERED BURGERS EQUATION *. R - 100,000 







